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STOCHASTIC GALERKIN METHODS FOR THE STEADY-STATE 
NAVIER-STOKES EQUATIONS * 

BEDRICH SOUSEDi'Kt AND HOWARD C. ELMAN* 


Abstract. We study the steady-state Navier-Stokes equations in the context of stochastic finite 
element discretizations. Specifically, we assume that the viscosity is a random field given in the form 
of a generalized polynomial chaos expansion. For the resulting stochastic problem, we formulate 
the model and linearization schemes using Picard and Newton iterations in the framework of the 
stochastic Galerkin method, and we explore properties of the resulting stochastic solutions. We 
also propose a preconditioner for solving the linear systems of equations arising at each step of 
the stochastic (Galerkin) nonlinear iteration and demonstrate its effectiveness for solving a set of 
benchmark problems. 


1. Introduction. Models of mathematical physics are typically based on partial 
differential equations (PDEs) that use parameters as input data. In many situations, 
the values of parameters are not known precisely and are modeled as random fields, 
giving rise to stochastic partial differential equations. In this study we focus on mod¬ 
els from fluid dynamics, in particular the stochastic Stokes and the Navier-Stokes 
equations. We consider the viscosity as a random field modeled as colored noise, and 
we use numerical methods based on spectral methods, specifically, the generalized 
polynomial chaos (gPC) framework [10, 14, 26, 27]. That is, the viscosity is given by 
a gPC expansion, and we seek gPC expansions of the velocity and pressure solutions. 

There is a number of reasons to motivate our interest in Navier-Stokes equations 
with stochastic viscosity. For example, the exact value of viscosity may not be known, 
due to measurement error, the presence of contaminants with uncertain concentra¬ 
tions, or of multiple phases with uncertain ratios. Alternatively, the fluid properties 
might be influenced by an external field, with applications for example in magnetohy¬ 
drodynamics. Specifically, we assume that the viscosity v depends on a set of random 
variables This means that the Reynolds number. 


Re(C) 


UL 

WY 


where U is the characteristic velocity and L is the characteristic length, is also stochas¬ 
tic. Consequently, the solution variables are random fields, and different realizations 
of the viscosity give rise to realizations of the velocities and pressures. As observed 
in [19], there are other possible formulations and interpretations of fluid flow with 
stochastic Reynolds number for example, where the velocity is fixed but the volume 
of fluid moving into a channel is uncertain so the uncertainty derives from the Dirichlet 
inflow boundary condition. 

We consider models of steady-state stochastic motion of an incompressible fluid 
moving in a domain D C Extension to three-dimensional models is straight¬ 
forward. We formulate the stochastic Stokes and Navier-Stokes equations using the 
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stochastic finite element method, assuming that the viscosity has a general probabil¬ 
ity distribution parametrized by a gPC expansion. We describe linearization schemes 
based on Picard and Newton iteration for the stochastic Galerkin method, and we 
explore properties of the solutions obtained, including a comparison of the stochas¬ 
tic Galerkin solutions with those obtained using other approaches, such as Monte 
Carlo and stochastic collocation methods [26]. Finally, we propose efficient hierar¬ 
chical preconditioners for iterative solution of the linear systems solved at each step 
of the nonlinear iteration in the context of the stochastic Galerkin method. Our ap¬ 
proach is related to recent work by Powell and Silvester [19]. However, besides using 
a general parametrization of the viscosity, our formulation of the stochastic Galerkin 
system allows straightforward application of state-of-the-art deterministic precondi¬ 
tioners by wrapping them in the hierarchical preconditioner developed in [22]. For 
alternative preconditioners see, e.g., [2, 11, 17, 18, 20, 23, 25]. Finally, we note that 
there exist related approaches based on stochastic perturbation methods [13], impor¬ 
tant developments also include reduced-order models such as [4, 24] , and an overview 
of existing methods for stochastic computational fluid dynamics can be found in the 
monograph [14]. 

The paper is organized as follows. In Section 2, we recall the deterministic steady- 
state Navier-Stokes equations and their discrete form. In Section 3, we formulate 
the model with stochastic viscosity, derive linearization schemes for the stochastic 
Galerkin formulation of the model, and explore properties of the resulting solutions 
for a set of benchmark problems that model the flow over an obstacle. In Section 4 we 
introduce a preconditioner for the stochastic Galerkin linear systems solved at each 
step of the nonlinear iteration, and in Section 5 we summarize our work. 

2. Deterministic Navier-Stokes equations. We begin by defining the model 
and notation, following [6]. For the deterministic Navier-Stokes equations, we wish to 
find velocity u and pressure p such that 

—vV^u + {u-V)u + Vp = f, (2.1) 

V-u = 0, (2.2) 

in a spatial domain D, satisfying boundary conditions 

u = g, on FDir, (2.3) 

v'S/u ■ fi — pfi = Q, on FNeui (2-4) 

where dD = Foir U Fnbu, and assuming sufficient regularity of the data. Dropping 
the convective term {u ■ V) u from (2.1) yields the Stokes problem 

Vp =/, (2.5) 

V-u = 0. (2.6) 

The mixed variational formulation of (2.1)-(2.2) is to find {u,p) G (Ve, Qd) such that 

V / Vu : VzT-l- / {u-Vu)v— / p{V-v)= / f-v, WvGVd, (2-7) 
J D J D J D J D 

f qiV-u)=0, yqGQo, (2.8) 

JD 
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where {Vd,Qd) is a pair of spaces satisfying the inf-sup condition and Ve is an 
extension of Ve containing velocity vectors that satisfy the Dirichlet boundary con¬ 
ditions [3, 6, 12]. 

Let c{z;u,v) = (z ■ Vm) • v. Because the problem (2.7)-(2.8) is nonlinear, it is 
solved using a linearization scheme in the form of Newton or Picard iteration, derived 
as follows. Consider the solution (u,p) of (2.7)-(2.8) to be given as u = + 6u^ 

and p = + 5p'^. Substituting into (2.7)-(2.8) and neglecting the quadratic term 

c{S'ir'-,Su^,v) gives 

ly [ : Vu-bc((5zr;zi”,u)-hc(?7”;(5zi”,u) - / (V • u) = i?” (u), 

Jd Jd 

f q(y.Sir)=r-{q), 

Jd 

where 

(v) = [ f-v-E [ ViT-.Vv- c{ir-, ir,v)+ [ p" (V • u), 

J D J D J D 

r^{q) = - [ qiV-iT). 

Jd 

Step n of the Newton iteration obtains (5tt",i5p") from (2.9)-(2.10) and updates the 
solution as 


( 2 . 11 ) 

( 2 . 12 ) 


(2.9) 

( 2 . 10 ) 


(2.13) 

(2.14) 

Step n of the Picard iteration omits the term c{5vP]vP,v) in (2.9), giving 


pU+l = pn + 


V [ N5Jr : Vv + c{Jr]5iN,v) 

Jd 


[ 6p^ (V • u) = i?" (v ), 
Jd 

[ qiV-Sir)= r" (g) . 
J D 


(2.15) 

(2.16) 


Consider the discretization of (2.1)-(2.2) by a div-stable mixed finite element 
method; for experiments discussed below, we used Taylor-Hood elements [6]. Let 
the bases for the velocity and pressure spaces be denoted and {piltu 

spectively. In matrix terminology, each nonlinear iteration entails solving a linear 
system 


F" 

■ Ju" ■ 


R" 

B 0 

Jp" 


j,n 


followed by an update of the solution 

u”+i = u”-H (5u", 

p"+i =p" + (5p". 


(2.17) 


(2.18) 

(2.19) 


For Newton’s method, F” is the Jacobian matrix, a sum of the vector-Laplacian 
matrix A, the vector-convection matrix N", and the Newton derivative matrix W", 


F” = A-^N" -hW”, 
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( 2 . 20 ) 








where 


^ — [^ab] : 

^ab — 

V [ Wcjib ■ V0a, 



Jd 

N" = K,] , 

II 

e 




Jd 

w” = [w:,] , 

w2b = 

[ ■ VW") • Jpa- 

Jd 


For Picard iteration, the Newton derivative matrix W" is dropped, and F” = A+N". 
The divergence matrix B is defined as 


— {^cd\ 5 ^cd — / * V^c) ■ 

JD 

The residuals at step n of both nonlinear iterations are computed as 


R’" 


■ f ■ 


pn 



■ u” ■ 

j.n 


. s 


B 

0 


p” 


( 2 . 21 ) 


( 2 . 22 ) 


where P" = A + N" and f is a discrete version of the forcing function of (2.1).^ 

3. The Navier-Stokes equations with stochastic viscosity. Let {Q.,F,V) 
represent a complete probability space, where 12 is the sample space, is a cr-algebra 
on 12 and P is a probability measure. We will assume that the randomness in the 
model is induced by a vector of independent, identically distributed (i.i.d.) random 
variables ^ such that ^ : 12 —F C K^. Let F denote the 

cr-algebra generated by and let /i(^) denote the joint probability density measure 
for The expected value of the product of random variables u and v that depend 
on ^ determines a Hilbert space Tp = {n,F^,fj,) with inner product 


{u,v) =E[uv] = J u {^) V dfi , (3.1) 

where the symbol E denotes mathematical expectation. 

3.1. The stochastic Galerkin formulation. The counterpart of the varia¬ 
tional formulation (2.7)-(2.8) consists of performing a Galerkin projection on the 
space Tp using mathematical expectation in the sense of (3.1). That is, we seek the 
velocity u, a random field in Ve ®Ty, and the pressure p £ Qd ®Tp, such that 


E 


/ i'Wu:\/v+ / (m-Vm)F— j p{y-v) =E f 
Jd Jd Jd J Ud 


VF S Ve ® Tp, 


E 


[ <7(V-F) 
Jd 


= 0 Vq £ Qd^Tt- 


The stochastic counterpart of the Newton iteration (2.9)-(2.10) is 


E 


[ lyVSiF :Vv + c{ir]Sir,v) + c{Sir-ir,v)- [ (5p”(V-F) 
J D J D 


E 


q{V -SJT) 


ID 


= i?" 


= r 


(3.2) 

(3.3) 


(3.4) 

(3.5) 


^Throughout this study, we use the convention that the right-hand sides of discrete systems 
incorporate Dirichlet boundary data for velocities. 
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where 


i?" (v) = E 
r” {q) = -E 


f-v— / ^ V?r* : ViT—c(?r*;?2",'i7) + / p" (V • zT) 


/£) 


g(V -zT) 


ID 


ID 


The analogue for Picard iteration omits c{5u^;u^,v) from (3.4): 


(3.6) 

(3.7) 


E 


[ uVSiT ■.\7v + c{ir;dir,v)- [ 6p^ {\7 ■ v) 
Jd J d 


= R^ 


(3.8) 


In computations, we will use a finite-dimensional subspace Tp C Tr spanned by a 
set of polynomials {V'z (■?)} that are orthogonal with respect to the density function /z, 
that is (z/’fcjz/’f) = Ski- This is referred to as the gPC basis; see [10, 27] for details 
and discussion. For Tp, we will use the space spanned by multivariate polynomials 
in of total degree P, which has dimension M = We will also 

assume that the viscosity is given by a gPC expansion 


M^-l 

^ t'z (a;) zAz (0 : (3-9) 

where {vi (x)} is a set of given deterministic spatial functions. 

3.2. Stochastic Galerkin finite element formnlation. We discretize (3.4) 
(or (3.8)) and (3.5) using div-stable finite elements as in Section 2 together with the 
gPC basis for Tp. For simplicity, we assume that the right-hand side / (x) and the 
Dirichlet boundary conditions (2.3) are deterministic. This means in particular that, 
as in the deterministic case, the boundary conditions can be incorporated into right- 
hand side vectors (specified as y below). Thus, we seek a discrete approximation of 
the solution of the form 

M-l M-1 

= EE (^) (0 ^ ^ (“C); (3.10) 

k—0 i—1 k—0 

M-l Np M-l 

p{x,0 ~ (3-11) 

k—0 j — 1 k—0 

The structure of the discrete operators depends on the ordering of the unknown 
coefficients {uik}, {Pjk}- We will group velocity-pressure pairs for each k, the index of 
stochastic basis functions (and order equations in the same way), giving the ordered 
list of coefficients 


Ul:N^,0,Pl:Np,0, Ul:Np,l,Pl-.Np,l, ■■■, Ul'.Np ,M-1, Pl:Np ,M-1 ■ (3-12) 

To describe the discrete structure, we first consider the stochastic version of the Stokes 
problem (2.5)-(2.6), where the convection term c(-; •, •) is not present in (3.4) and (3.8). 
The discrete stochastic Stokes operator is built from the discrete components of the 
vector-Laplacian 


A£= [ae^ab], at.ab = { (^) 


'D 


5 


^ = 1, (3.13) 








which are incorporated into the block matrices 


5o = 


Ao 

B 0 


5^ 


Af 0 
0 0 ’ 


- 1 . 


(3.14) 


These operators will be coupled with matrices arising from terms in Tp, 

iii = [hijk\, = E , € = 0,... - 1, j,k = 0,...,M -1. 

(3.15) 

Combining the expressions from (3.13), (3.14) and (3.15) and using the ordering (3.12) 
gives the discrete stochastic Stokes system 


^ ^ Hf (8)5^^ V = y, (3.16) 

where ® corresponds to the matrix Kronecker product. The unknown vector v corre¬ 
sponds to the ordered list of coefficients in (3.12) and the right-hand side is ordered 
in an analogous way. Note that Hq is the identity matrix of order M. 

Remark 3.1. With this ordering, the coefficient matrix contains a set of M 
block 2x2 matrices of saddle-point structure along its block diagonal, given by 


M„-l 

So -h hijjS^, j = 0,..., M — 1. 
e=i 


This enables the use of existing deterministic solvers for the individual diagonal blocks. 
An alternative ordering based on the blocking of all velocity coefficients followed by 
all pressure coefficients, considered in [19], produces a matrix of global saddle-point 
structure. 

The matrices arising from the linearized stochastic Navier-Stokes equations aug¬ 
ment the Stokes systems with stochastic variants of the vector-convection matrix and 
Newton derivative matrix appearing in (2.20). In particular, at step n of the nonlinear 
iteration, let (x) be the £th term of the velocity iterate (as in the expression on the 
right in (3.10) for k = £), and let 


N? = Kab] , 

n7,ab= [ (z^ • V4) • </>a, 

Jd 


wr = K,,], 

Wlab= f {c[b-Vffil)-<[a. 

JD 


Then the analogues of (3.13)-(3.14) 

are 


F'f = Ae + 'N'f 

F'f = Ae + 'N^, 

for the stochastic Newton method 

for stochastic Picard iteration. 

(3.17) 

(3.18) 


so for Newton’s method 


■pn _ 
■^0 — 


FO B"^ 

B 0 


pri _ 
•J 9 — 


F^ 0 
0 0 


(3.19) 


and as above, for Picard iteration the Newton derivative matrices {W"} are dropped. 
Note that £ = 0,..., M — 1 here, where M = max (M, M,^). (In particular, if > M, 
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we set = W" = 0 for £ = M + 1,..., — 1.) Step n of the stochastic nonlinear 

iteration entails solving a linear system and updating, 


M-l 




= 7^", 


^n+l = v" +(5v", 


(3.20) 


where 


= y- 


M-l 

E 


t=0 


(3.21) 


v" and (5v" are vectors of current velocity and pressure coefficients and updates, 
respectively, ordered as in (3.12), y is the similarly ordered right-hand side determined 
from the forcing function and Dirichlet boundary data, and 

Af + 0 \ 


'Po = 


-NS 


B 


0 


note that the (1, l)-blocks here are as in (3.18). 

3.3. Sampling methods. In experiments described below, we compare some 
results obtained using stochastic Galerkin methods to those obtained from Monte 
Carlo and stochastic collocation. We briefly describe these approaches here. 

Both Monte Carlo and stochastic collocation methods are based on sampling. 
This entails the solution of a number of mutually independent deterministic problems 
at a set of sample points which give realizations of the viscosity (3.9). That 

is, a realization of viscosity Jz gives rise to deterministic functions 

and on D that satisfy the standard deterministic Navier-Stokes equations, 

and to finite-element approximations {x ), 

In the Monte Carlo method, the Nmc sample points are generated randomly, 
following the distribution of the random variables and moments of the solution 
are obtained from ensemble averaging. For stochastic collocation, the sample points 
consist of a set of predetermined collocation points. This approach derives from a 
methodology for performing quadrature or interpolation in multidimensional space 
using a small number of points, a so-called sparse grid [8, 16]. There are two ways to 
implement stochastic collocation to obtain the coefficients in (3.10)-(3.11), either by 
constructing a Lagrange interpolating polynomial, or, in the so-called pseudospectral 
approach, by performing a discrete projection into Tp [26]. We use the second ap¬ 
proach because it facilitates a direct comparison with the stochastic Calerkin method. 
In particular, the coefficients are determined using a quadrature 

q=l q=l 

where are collocation (quadrature) points, and are quadrature weights. We 
refer, e.g., to [14] for an overview and discussion of integration rules. 

3.4. Example: flow around an obstacle. In this section, we present results 
of numerical experiments for a model problem given by a flow around an obstacle 
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Fig. 3.1. Finite element mesh for the flow around an obstacle problem. 


in a channel of length 12 and height 2. The viscosity (3.9) was taken to be a trun¬ 
cated lognormal process with mean values vq = 1/50 or 1/150, which corresponds 
to mean Reynolds numbers Reo = 100 or 300, respectively, and its representation 
was computed from an underlying Gaussian random process using the transformation 
described in [9]. That is, for £ = 0,..., M^, — 1, ij)i (^) is the product of N univariate 
Hermite polynomials, and denoting the coefficients of the Karhunen-Loeve expansion 
of the Gaussian process by gj (x) and rjj = — gj, j = 1,..., N, the coefficients in 

the expansion (3.9) are computed as 


(x) = E [ipe {r])] exp 


9o (a:) + 


I i9j {x) f 


The covariance function of the Gaussian held, for points Xi = {xi,yi), X 2 = (x 2 , 92 ) G 
H, was chosen to be 


C{Xi,X 2 ) = o-|exp 


\X2 - Xi\ 


L 


X 


\y2 - yih 

Ly J’ 


where and Ly are the correlation lengths of the random variables i = 1,..., X, 
in the x and y directions, respectively, and ag is the standard deviation of the Gaussian 
random held. The correlation lengths were set to be equal to 25% of the width and 
height of the domain, i.e. = 3 and Ly — 0.5. The coefficient of variation of the 
lognormal held, dehned as CoV = OvIvq where Uv is the standard deviation, was 10% 
or 30%. The stochastic dimension was N = 2. The lognormal formulation was chosen 
to enable exploration of a general random held in which the viscosity guaranteed to be 
positive. See [1] for an example of the use of this formulation for dihusion problems 
and [28] for its use in models of porous media. 

We implemented the methods in Matlab using IFISS 3.3 [5]. The spatial dis¬ 
cretization uses a stretched grid, discretized by 1520 Taylor-Hood hnite elements; the 
domain and grid are shown in Figure 3.1. There are 12,640 velocity and 1640 pressure 
degrees of freedom, the degree used for the polynomial expansion of the solution was 
P = 3, and the degree used for the expansion of the lognormal process was 2P, which 
ensures a complete representation of the process in the discrete problem [15]. With 
these settings, M = 10 and = M = 28, and is of order 10 in (3.20). 

Gonsider hrst the case of Reo = 100 and CoV — 10%. Figure 3.2 shows the mean 
horizontal and vertical components of the velocity and the mean pressure (top), and 
the variances of the same quantities (bottom). It can be seen that there is symmetry 
in all the quantities, the mean values are essentially the same as we would expect 
in the deterministic case, and the variance of the horizontal velocity component is 
concentrated in two “eddies” and is larger than the variance of the vertical velocity 
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Fig. 3.2. Mean horizontal and vertical velocities and pressure (top) and variances of the same 
quantities (bottom), for Re = 100 and CoV = 10%. 
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component. Figure 3.3 illustrates values of several coefficients of expansion (3.10) of 
the horizontal velocity. All the coefficients are symmetric, and as the index increases 
they become more oscillatory and their values decay. We found the same trends for the 
coefficients of the vertical velocity component and of the pressure. Our observations 
are qualitatively consistent with numerical experiments of Powell and Silvester [19]. 



Fig. 3.3. Coefficients 1 — 4 o/ the gPC expansion of the horizontal velocity, Reo = 100 and 
CoV = 10%. 


We also tested (the same) Reo = 100 with increased CoV = 30%. We found that 
the mean values are essentially the same as in the previous case; Figure 3.4 shows 
the variances, which display the same qualitative behavior but have values that are 
approximately 10 times larger than for the case CoV = 10%. 

A different perspective on the results is given in Figure 3.5, which shows estimates 
of the probability density function (pdf) for the horizontal velocity at two points in 
the domain, (4.0100, —0.4339) and (4.0100,0.4339). These are locations at which large 
variances of the solution were seen, see Figures 3.2 and 3.4. The results were obtained 
using Matlab’s ksdensity function. It can be seen that with the larger value of CoV, 
the support of the velocity pdf is wider, and except for the peak values, for fixed CoV 
the shapes of the pdfs at the two points are similar, indicating a possible symmetry 
of the stochastic solution. For this benchmark, we also obtained analogous data using 
the Monte Carlo and collocation sampling methods; it can be seen from the figure 
that these methods produced similar results. ^ 

Next, we consider a larger value of the mean Reynolds number, Reo = 300. 
Figure 3.6 shows the means and variances for the velocities and pressure for CoV = 
10%. It is evident that increased Reo results in increased values of the mean quantities, 


^The results for Monte Carlo were obtained using 10® samples, and those for collocation were 
found using a Smolyak sparse grid with Gauss-Hermite quadrature and grid level riq = 4. 
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Fig. 3.4. Variances of velocity components and pressure for Reo = 100 and CoV = 30%. 


but they are again similar to what would be expected in the deterministic case. The 
variances exhibit wider eddies than for Rep = 100, and in this case there is only one 
region of the largest variance in the horizontal velocity, located just to the right of 
the obstacle; this is also a region with increased variance of the pressure. 

In similar tests for the larger value CoV = 30%, we found that the mean values are 
essentially the same as for CoV = 10%, and Figure 3.7 shows the variances of velocities 
and pressures. From the figure it can be seen that the variances are qualitatively the 
same but approximately 10 times larger than for CoV = 10%, results similar to those 
found for Reg = 100. 

As above, we also examined estimated probability density functions for the ve¬ 
locities and pressures at a specified point in the domain, in this case, at the point 
(3.6436,0), taken again from the region in which the solution has large variance. Fig¬ 
ure 3.8 shows these pdf estimates, from which it can be seen that the three methods 
for handling uncertainty are in close agreement for each of the two values of CoV. 

Finally, we show in Figure 3.9 the results of one additional experiment with 
Reg = 300 and CoV = 10%, where estimated pdfs of the first velocity component Ux 
were computed at three points near the inflow boundary, {x,y) = (0.5,0), (1,0) and 
(1.5,0). These plots show some effects of spatial accuracy. The results for all methods 
were identical, and so only one representative pdf is shown. The images on the top 
and bottom left show results for a uniform mesh of width /i = 1/8 and for two refined 
meshes in which the horizontal mesh width to the left of the obstacle is reduced to 
h/2 and h/A. The image in the bottom right provides a more detailed view of the 
fine-grid results; in this image, the width of the horizontal window is the same (0.01) 
for the three subplots but the vertical heights are different. Several trends are evident: 

• The pdfs at points nearer the boundary are much narrower. 

• When the spatial accuracy is low, the methods miss some features of the pdf. 
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Fig. 3.5. Estimated pdfs of the velocities with Reo = 100, CoV = 10% (left) and 30% (right) 
at the points with coordinates (4.0100,-0.4339) (top) and (4.0100,0.4339) (bottom). 


In particular the methods produce a pdf that captures its “spiky” nature near 
the inflow, but the location of the spike is not correct. 

We believe these effects stem from the fact that the inflow boundary is deterministic 
(with My = 1 at y = 0), and the effects of variability in the viscosity are felt less 
strongly at points near the inflow boundary. At points further from the deterministic 
inflow boundary, these effects and differences become less dramatic. 

3.5. Nonlinear solvers. We briefly comment on the nonlinear solution algo¬ 
rithm used to generate the results of the previous section. The nonlinear solver was 
implemented by modifying the analogue for deterministic systems in IFISS. It uses 
a hybrid strategy in which an initial approximation is obtained from solution of the 
stochastic Stokes problem (3.16), after which several steps of Picard iteration (equa¬ 
tion (3.20) with (Fg specified using (3.19) and (3.18)) are used to improve the solution, 
followed by Newton iteration {Ft, from (3.17)). A convergent iteration stopped when 
the Euclidian norm of the algebraic residual (3.21) satisfied |17^"||2 < ellj/lb where 
£ = 10“® and y is as in (3.16). 

In the experiments described in Section 3.4, we used values of the Reynolds num¬ 
ber, Re = 100 and 300, and for each of these, two values of the coefficient of variation, 
CoV = 10% and 30%. We list here the numbers of steps leading to convergence of 
the nonlinear algorithms that were used to generate the solutions discussed above. 
Direct solvers were used for the linear systems; we discuss a preconditioned iterative 
algorithm in Section 4 below. 

Re = 100, CoV = 10%: 6 Picard steps 1 Newton step(s) 

Re = 100, CoV = 30%: 6 2 
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Fig. 3.6. Mean horizontal and vertical velocities and pressure (top) and variances of the same 
quantities (bottom), for Reo = 300 and CoV = 10%. 
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Fig. 3.7. Variances of velocity components and pressure for Reo = 300 and CoV = 30%. 


Re = 300, CoV = 10%: 20 1 

Re = 300, CoV = 30%: 20 2 

Thus, a larger CoV (larger standard deviation of the random field determining uncer¬ 
tainty in the process) leads to somewhat larger computational costs. For Re = 300, 
the nonlinear iteration was not robust with 6 initial Picard steps (for the stochastic 
Galerkin method as well as the sampling methods); 20 steps was sufficient. 

We also explored an inexact variant of these methods, in which the coefficient 
matrix of (3.20) for the Picard iteration was replaced by the block diagonal matrix 
Ho ® obtained from the mean coefficient. For CoV = 10%, with the same number 
of (now inexact) Picard steps as above (6 for Re = 100 and 20 for Re = 300), this 
led to just one extra (exact) Newton step for Re = 100 and no additional steps for 
Re = 300. On the other hand, for CoV = 30%, this inexact method failed to converge. 

4. Preconditioner for the linearized systems. The solution of the linear 
systems required during the course of the nonlinear iteration is a computationally 
intensive task, and use of direct solvers may be prohibitive for large problems. In this 
section, we present a preconditioning strategy for use with Krylov subspace methods 
to solve these systems, and we compare its performance with that of several other 
techniques. The new method is a variant of the hierarchical Gauss-Seidel precondi¬ 
tioner developed in [22], 

4.1. Structure of the matrices and the preconditioner. We first recall the 
structure of the matrices {H,} of (3.15). More comprehensive overviews of these 
matrices can be found in [7, 15]. The matrix structure can be understood through 
knowledge of the coefficient matrix cp = where j,k = 0,..., M — 1. 

The block sparsity structure depends on the type of coefficient expansion in (3.9). If 
only linear terms are included, that is tpi = i = l,---,.^, then the coefficients 
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Fig. 3.8. Estimated pdfs of the velocities (top), Uy (middle), and pressurep (bottom) with 
Reo = 300 and CoV = 10% (left) and 30% (right) at the point with coordinates (3.6436,0). 


hijk = E yield a Galerkin matrix with a block sparse structure. In the more 

general case, = E [tjtiipjipk] and the stochastic Galerkin matrix becomes fully 

block dense. In either case, for fixed £ and a set of degree (P polynomial expansions, 
with 1 < y < P, the corresponding coefficient matrix cy has a hierarchical structure 


c-p 


cp-i 

bp dp 




Now, let Ap denote the global stochastic Galerkin matrix corresponding to either a 
Stokes problem (3.16) or a linearized system (3.20); we will focus on the latter system 
in the discussion below. The matrix Ap also has a hierarchical structure 


Ap = 


Ap-i Bp 
Cp Vp ’ 


T = 1,...,P, 
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(4.1) 


















































Fig. 3.9. Estimated pdfs of the velocity at points (0.5,0) (top left) (1,0) (top right) and 
(1.5,0) (bottom left), with Reo = 300 and CoV = 10% and three meshes. Bottom right: detailed 
depiction of the pdfs for each grid point obtained from the finest mesh. 


where Ao is the matrix of the mean, derived from vq in (3.9). This hierarchical 
structure is shown in the left side of Figure 4.1. 

We will write vectors with respect to this hierarchy as 


a;(o:y) 


a;(o) 

Xil) 


. . 


where includes all indices corresponding to polynomial degree q, blocked by spatial 
ordering determined by (3.12). With this notation, the global stochastic Galerkin 
linear system has the form 


ApX(^0,p') — /(0:P). 


(4.2) 


To formulate the preconditioner for (4.2), we let Ao represent an approximation 
of Ao and Py represent an approximation of Py. In particular, let 


Vy = 


Ao 


Ao _ 


(4.3) 
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where the nunAer of diagonal blocks is given by IP. We will need the action of 
the inverse of or an approximation to it, which can be obtained using an LU- 
factorization of Ao, or using some preconditioner for Aq., or using a Krylov subspace 
solver. A preconditioner V : wi^o-.y) U(o:a>) for (4.2) is then defined as follows: 

Algorithm 4.1. [Approximate hierarchical Gauss-Seidel preconditioner (ahGS)] 
Solve (or solve approximately) 

AoT(o)=W(o), (4.4) 

and, for T= 1,... P, solve (or solve approximately) 

i>yV(y) = (w(y) - Cg>U( 0 :T-i)) . (4.5) 

The cost of preconditioning can be reduced further by truncating the matrix- 
vector (MATVEC) operations used for the multiplications by the submatrices Cy 
in (4.5). The idea is as follows. The system (4.2) can be written as 

M-l M„-l 

5^ hejkJ^iXj = fk, k = 0,...,M -1, (4.6) 

j-o e=o 

and the MATVEC with Ap is given by 


M-l M„-l 

Tj — ^ ^ ^ ^ (^■^) 

k=0 £=0 


where the indices j, k G — 1} correspond to nonzero blocks in Ap. The 

truncated MATVEC is an inexact evaluation of (4.7) proposed in [22, Algorithm 1], 
in which the summation over £ = 0,... .M^ — 1 is replaced by summation over a subset 
Ait C {0,..., Ml, — 1}. Figure 4.1 shows the hierarchical structure of the matrix and 
of the ahGS preconditioning operator. Both images in the figure correspond to the 
choice P = 3, so that the hierarchical preconditioning operation (4.4)-(4.5) requires 


four steps. Because N = 4, the matrix block size is M = 


JV-I-F 

P 


= 35. 


The 


block-lower-triangular component of the image on the right in Figure 4.1 shows the 
hierarchical structure of the ahGS preconditioning operator with truncation. For the 


matrix in the left panel, M,, 


N + 2P 
2P 


= 210, but the index set Mt includes 


terms with indices at most M — 1 in the accumulation of sums used for C-p. These two 
heuristics, approximation by (4.3) and the truncation of MATVECs, significantly im¬ 
prove the sparsity structure of the preconditioner, on both the block diagonal (through 
the first technique) and the block lower triangle (through the second). In the next 
section, we will also consider truncated MATVEGs with smaller maximal indices. 


4.2. Numerical experiments. In this section, we describe the results of exper¬ 
iments in which the ahGS preconditioner is used with GMRES to solve the systems 
arising from both Picard and Newton iteration.^ Unless otherwise specified, in the 


®The preconditioning operator may be nonlinear, for example, if the block solves in (4.4)—(4.5) 
are performed approximately using Krylov subspace methods, so that care must be made in the 
choice of the Krylov subspace iterative method. For this reason, we used a flexible variant of the 
GMRES method [21]. 
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Fig. 4.1. Hierarchical structure of the stochastic Galerkin matrix (4-1) (left) and splitting 
operator L-\-diag{D) for the approximate hierarchical Gauss-Seidel preconditioner (ahGS) with the 
truncation of the MATVEC (right). 


experiments discussed below, the block-diagonal matrices of (4.3), (4.4), (4.5) are 
taken to be Aq = Aq, which corresponds to the first summand Hg (8 ).Fq in (3.20) and 
represents an approximation to the diagonal blocks of Vy; see also Remark 3.1. We 
also compare the performance of ahGS preconditioning with several alternatives: 

• Mean-based preconditioner (MB) [17, 18], where the preconditioning operator 
is the block-diagonal matrix Hg G Tq = J G -Fq derived from the mean of i'. 

• The Kronecker-product preconditioner [25] (denoted K), given by Hg ® Fq, 
where Hg is chosen to to minimize a measure of the difference Ap — IIo^Fq. 

• Block Gauss-Seidel preconditioner (bGS), in which the preconditioning oper¬ 
ation entails applying (Py -|- L)~^ determined using the block lower triangle 
of A-p, that is with the approximation of the diagonal blocks by Aq but 
without MATVEG truncation. This is an expensive operator but enables an 
understanding of the impact of various efforts to make the preconditioning 
operator more sparse. 

• ahGS(PGD), a modification of ahGS in which the block-diagonal matrices of 
(4.4)-(4.5) are replaced by the pressure conveetion-diffusion (PCD) approxi¬ 
mation [6] to the mean matrix Fq ■ 

• ahGS(PCD-it), in which the block-diagonal solve of (4.5) is replaced by an 
approximate solve determined by twenty steps of a PCD-preconditioned GM- 
RES iteration. For this strategy, Ao in (4.4) and the blocks of T>-j> in (4.5) 
are taken to be the complete sum of (3.20), but the approximate solutions to 
these systems are obtained using a preconditioned inner iteration. 

Four of the strategies, the ahGS, mean-based, Kronecker-product and bGS pre¬ 
conditioners, require the solution of a set of block-diagonal systems with the structure 
of a linearized Navier-Stokes operator of the form given in (2.17). We used direct 
methods for these computations. All results presented are for Re = 100; performance 
for Re = 300 were essentially the same. We believe this is because of the exact solves 
performed for the mean operators. 

The results for the first step of Picard iteration are in Tables 4.1-4.4. All tests 
started with a zero initial iterate and stoppped when the residual r^^'> = /(g:P) — 

•^pfo-p) iterate satisfied llr(^^j |2 < 10“®]]l/(g:P)II 2 in the Euclidian 
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norm. With other parameters fixed and no truncation of the MATVEC, Table 4.1 
shows the dependence of GMRES iterations on the stochastic dimension N, Table 4.2 
shows the dependence on the degree of polynomial expansion P, and Table 4.3 shows 
the dependence on the coefficient of variation CoV. It can be seen that the num¬ 
bers of iterations with the ahGS preconditioner are essentially the same as with the 
block Gauss-Seidel (bGS) preconditioner, and they are much smaller compared to 
the mean-based (MB) and the Kronecker product preconditioners. When the ex¬ 
act solves with the mean matrix are replaced by the mean-based modified pressure- 
convection-diffusion (PCD) preconditioner for the diagonal block solves, the iterations 
grow rapidly. On the other hand, if the PCD preconditioner is used as part of an in¬ 
ner iteration as in ahGS(PCD-it), then good performance is recovered. This indicates 
that a good preconditioner for the mean matrix is an essential component of the global 
preconditioner for the stochastic Galerkin matrix. 

Table 4.4 shows the iteration counts when the MATVEC operation is truncated 
in the action of the preconditioner. Truncation decreases the cost per iteration of the 
computation, and it can be also seen that performance can actually be improved. Eor 
example, with it = 2, the number of iterations is the smallest. Moreover there are 
only 21 nonzeros in the lower triangular part of the sum of coefficient matrices {H^} 
(each of which has order 10) used in the MATVEC with if < 2, compared to 63 
nonzeros when no truncation is used; there are 203 nonzeros in the set of 28 full 
matrices {H^}. 


Table 4.1 

For Picard step: dependence on stochastic dimension N of GMRES iteration counts, for various 
preconditioners, with polynomial degree P = S and coefficient of variation CoV = 30%. M is the 
block size of the stochastic Galerkin matrix, the number of terms in (3.9) and ngdof the size of 
the stochastic Galerkin matrix. 


N 

M 

M, 

ngdof 

MB 

K 

bCS 

ahCS 

ahGS(PCD) 

ahGS(PCD-it) 

1 

4 

7 

57,120 

63 

36 

30 

30 

201 

29 

2 

10 

28 

142,800 

102 

66 

59 

54 

357 

48 

3 

20 

84 

285,600 

145 

109 

88 

82 

553 

73 


Table 4.2 

For Picard step: dependence on polynomial degree P of GMRES iteration counts, for various 
preconditioners, with stochastic dimension N = 3 and coefficient of variation CoV = 30%. Other 
headings are as in Table 4G. 


P 

M 


ngdof 

MB 

K 

bCS 

ahCS 

ahGS(PCD) 

ahGS(PCD-it) 

1 

4 

10 

57,120 

26 

22 

12 

11 

144 

13 

2 

10 

35 

142,800 

63 

49 

29 

26 

300 

29 

3 

20 

84 

285,600 

145 

109 

88 

82 

553 

73 


Results for the first step of Newton iteration, which comes after six steps of 
Picard iteration, are summarized in Tables 4.5-4.8. As above, the first three tables 
show the dependence on stochastic dimension N (Table 4.5), polynomial degree P 
(Table 4.6), and coefficient of variation CoV (Table 4.7). It can be seen that all 
iteration counts are higher compared to corresponding results for the Picard iteration, 
but the other trends are very similar. In particular, the performances of the ahCS 
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Table 4.3 

For Picard step: dependence on coefficient of variation CoV of GMRES iteration counts, for 
various preconditioners, with stochastic dimension N = 2 and polynomial degree P = 3. Other 
headings are as in Table 4-1, 


CoV{Vo) 

MB 

K 

bGS 

ahGS 

ahGS(PCD) 

ahGS(PCD-it) 

10 

16 

14 

7 

7 

186 

10 

20 

36 

27 

17 

16 

259 

17 

30 

102 

66 

59 

54 

357 

48 


Table 4.4 

For Picard step: number of GMRES iterations when the preconditioners use truncated 
MATVEG. The matrices corresponding to higher degree expansion of the coefficient than It are 
dropped from the action of the preconditioner. Here N = 2, P = 3, Mt is the number of terms used 
in the inexact (truncated) evaluation of the MATVEG (4-^)! nnz(ct) is the number of nonzeros in 
the lower triangular parts of the coefficient matrices (3.15) after the truncation. Other headings are 
as in Table f-l- 


it 

0 

1 

2 

3 

6 

setup Mt 

1 

3 

6 

10 

28 

nnz{c£) 

0 

12 

21 

43 

63 

bGS 

10 

16 

9 

7 

7 

7 

CoV (%) 20 

36 

19 

15 

17 

17 

30 

102 

55 

43 

57 

58 

ahGS 

10 

16 

9 

7 

7 

7 

CoV (%) 20 

36 

19 

14 

16 

16 

30 

102 

55 

35 

55 

54 


and bGS preconditioners are comparable except for the case when N = S, and P = 3 
(the last rows in Tables 4.5 and 4.6). Nevertheless, checking results in Table 4.8, 
which shows the effect of truncation, it can be seen that with the truncation of the 
MATVEG the iteration counts of the ahGS and bGS preconditioners can be further 
improved. Indeed, running one more experiment for the aforementioned case when 
V = 3, and P = 3, it turns out that with it = 2 the number of iterations with the 
ahGS and bGS preconditioners are 118 and 120, respectively and with it = i they are 
101 and 104, respectively. That is, the truncation leads to fewer iterations also in this 
case, and the performance of ahGS and bGS preconditioners is again comparable. 

Table 4.5 

For Newton step: dependence on stochastic dimension N of GMRES iteration counts, for 
various preconditioners, with polynomial degree P = 3 and coefficient of variation CoV = 30%. 
Other headings are as in Table 4-1. 


N 

M 

M, 

ngdof 

MB 

K 

bGS 

ahGS 

ahGS(PCD) 

ahGS(PCD-it) 

1 

4 

7 

57,120 

73 

42 

32 

32 

309 

40 

2 

10 

28 

142,800 

126 

80 

77 

69 

546 

63 

3 

20 

84 

285,600 

235 

170 

128 

151 

1011 

117 


Finally, we briefly discuss computational costs. For any preconditioner, each 
GMRES step entails a matrix-vector product by the coefficient matrix. For viscosity 
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Table 4.6 

For Newton step: dependence on polynomial degree P of GMRES iteration counts, for various 
preconditioners, with stochastic dimension N = 3 and coefficient of variation CoV = 30%. Other 
headings are as in Table f.l. 


p 

M 

M, 

ngdof 

MB 

K 

bGS 

ahGS 

ahGS(PCD) 

ahGS(PGD-it) 

1 

4 

10 

57,120 

27 

23 

12 

11 

219 

32 

2 

10 

35 

142,800 

68 

55 

33 

32 

450 

42 

3 

20 

84 

285,600 

235 

170 

128 

151 

1011 

117 


Table 4.7 

For Newton step: dependence on coefficient of variation CoV of GMRES iteration counts, for 
various preconditioners, with stochastic dimension N = 2 and polynomial degree P = 3. Other 
headings are as in Table f.l. 


CoV{%) 

MB 

K 

bGS 

ahGS 

ahGS(PGD) 

ahGS(PCD-it) 

10 

17 

15 

8 

8 

322 

28 

20 

40 

30 

20 

19 

379 

35 

30 

126 

80 

77 

69 

546 

63 


given by a general probability distribution, this will typically involve a block-dense 
matrix, and, ignoring any overhead associated with increasing the number of GMRES 
steps, this will be the dominant cost per step. The mean-based preconditioner requires 
the action of the inverse of the block-diagonal matrix / This has relatively small 

amount of overhead once the factors of are computed. The ahGS preconditioner 
without truncation effectively entails a matrix-vector product by the block lower- 
triangular part of the coefficient matrix, so its overhead is bounded by 50% of the cost 
of a multiplication by the coefficient matrix. This is an overestimate because it ignores 
the approximation of the block diagonal and the effect of truncation. For example, 
consider the case with stochastic dimension N = 2 and polynomial expansions of 
degree P = 3 for the solution and 2P = 6 for the viscosity; this gives M = 10 and 

= 28. In Tables 4.4 and 4.8, Mf indicates how many matrices are used in the 
MATVEG operations and nnz{ci) is the number of nonzeros in the sum of the lower 
triangular parts of the coefficient matrices {H^ | £ = 0,..., Mt — 1}. With complete 
truncation, ft = 0, and ahGS reduces to the mean-based preconditioner. With no 
truncation, ft = 6, and because the number of nonzeros in {H^} is 203, the overhead 
of ahGS is 63/203, less than 30% of the cost of multiplication by the coefficient 
matrix. If truncation is used, in particular when the iteration count is the lowest 
{ft = 2), the overhead is only 21/203 = 10%. Note that with increasing stochastic 
dimension and degree of polynomial expansion, the savings will be higher because 
the ratio of the sizes of the blocks C-y/Vy decreases as P increases, see (4.1). Last, 
the mean-based preconditioner is embarrassingly parallel; the ahGS preconditioner 
requires P-l-1 sequential steps, although each of these steps is also highly parallelizable. 
The Kronecker preconditioner is more difficult to assess because it does not have 
block-diagonal structure, and we do not discuss it here. 

5. Conclusion. We studied the Navier-Stokes equations with stochastic viscos¬ 
ity given in terms of polynomial chaos expansion. We formulated the stochastic 
Galerkin method and proposed its numerical solution using a stochastic versions of 
Picard and Newton iteration, and we also compared its performance in terms of ac- 
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Table 4.8 

For Newton step: number of GMRES iterations when the preconditioners use the truncation of 
the MATVEC. Headings are as in Table 


it 

0 

1 

2 

3 

6 

setup Mt 

1 

3 

6 

10 

28 

nnz{ci) 

0 

12 

21 

43 

63 

bGS 

10 

17 

10 

8 

8 

8 

CoV (%) 20 

40 

22 

18 

20 

20 

30 

126 

68 

57 

75 

77 

ahGS 

10 

17 

10 

8 

8 

8 

GoV (%) 20 

40 

22 

17 

19 

19 

30 

126 

68 

45 

70 

69 


curacy with that of stochastic collocation and Monte Carlo method. Finally, we 
presented a methodology of Gauss-Seidel hierarchical preconditioning with approxi¬ 
mation using the mean-based diagonal block solves and a truncation of the MATVEC 
operations. The advantage of this approach is that neither the matrix nor the precon¬ 
ditioner need to be formed explicitly, and the ingredients include only the matrices 
from the polynomial chaos expansion and a good preconditioner for the mean-value 
deterministic problem, it allows an obvious parallel implementation, and it can be 
written as a “wrapper” around existing deterministic code. 
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